home *** CD-ROM | disk | FTP | other *** search
/ Over 1,000 Windows 95 Programs / Over 1000 Windows 95 Programs (Microforum) (Disc 1).iso / 1249 / eigenval.t < prev    next >
Text File  |  1997-04-18  |  14KB  |  790 lines
  1. %
  2. % "eigenval.t" computes algebraic eigenvalues 
  3. % using QR decomposition
  4. %
  5. % adapted from: 
  6. % "C Tools for Scientists and Engineers" by L. Baker
  7. %
  8. %   Sample program for the T Interpreter by:
  9. %
  10. %   Stephen R. Schmitt
  11. %   962 Depot Road
  12. %   Boxborough, MA 01719
  13. %
  14.  
  15. var macheps : real := 0.00000001
  16. const DIM : int := 3
  17.  
  18. type rmatrix : array[DIM,DIM] of real
  19. type ivector : array[DIM] of int
  20. type rvector : array[DIM] of real
  21.  
  22. var A : rmatrix 
  23.  
  24. program
  25.  
  26.     var wr, wi : rvector
  27.     var z : rmatrix
  28.     var i : int
  29.     var iwork, cnt : ivector
  30.  
  31.     put "array 1:"
  32.  
  33.     A[0,0] := 2
  34.     A[0,1] := 2
  35.     A[0,2] := 1
  36.  
  37.     A[1,0] := 1
  38.     A[1,1] := 3
  39.     A[1,2] := 1
  40.  
  41.     A[2,0] := 1
  42.     A[2,1] := 2
  43.     A[2,2] := 2
  44.     print_matrix( A )
  45.  
  46.     eigenvv( A, iwork, z, wr, wi, cnt )
  47.  
  48.     put ""
  49.     put "eigenvalues for array 1:"
  50.     put ""
  51.  
  52.     for i := 0...DIM-1 do
  53.         put complex_str( wr[i], wi[i] ), "   it = ", cnt[i]
  54.     end for
  55.  
  56.  
  57.  
  58.     put ""
  59.     put ""
  60.     put "array 2:"
  61.  
  62.     A[0,0] := 2
  63.     A[0,1] := 4
  64.     A[0,2] := 4
  65.  
  66.     A[1,0] := 0
  67.     A[1,1] := 3
  68.     A[1,2] := 1
  69.  
  70.     A[2,0] := 0
  71.     A[2,1] := 1
  72.     A[2,2] := 3
  73.     print_matrix( A )
  74.  
  75.     eigenvv( A, iwork, z, wr, wi, cnt )
  76.  
  77.     put ""
  78.     put "eigenvalues for array 2:"
  79.     put ""
  80.  
  81.     for i := 0...DIM-1 do
  82.         put complex_str( wr[i], wi[i] ), "   it = ", cnt[i]
  83.     end for
  84.  
  85.  
  86.     
  87.     put ""
  88.     put ""
  89.     put "array 3:"
  90.  
  91.     A[0,0] :=  5
  92.     A[0,1] := -4
  93.     A[0,2] := -7
  94.  
  95.     A[1,0] := -4
  96.     A[1,1] :=  2
  97.     A[1,2] := -4
  98.  
  99.     A[2,0] := -7
  100.     A[2,1] := -4
  101.     A[2,2] :=  5
  102.     print_matrix( A )
  103.  
  104.     put ""
  105.     put "eigenvalues for array 3:"
  106.     put ""
  107.  
  108.     eigenvv( A, iwork, z, wr, wi, cnt )
  109.  
  110.     for i := 0...DIM-1 do
  111.         put complex_str( wr[i], wi[i] ), "   it = ", cnt[i]
  112.     end for
  113.  
  114.  
  115.     
  116.     put ""
  117.     put ""
  118.     put "array 4:"
  119.  
  120.     A[0,0] :=  1
  121.     A[0,1] := -1
  122.     A[0,2] := -1
  123.  
  124.     A[1,0] :=  1
  125.     A[1,1] := -1
  126.     A[1,2] :=  0
  127.  
  128.     A[2,0] :=  1
  129.     A[2,1] :=  0
  130.     A[2,2] := -1
  131.     print_matrix( A )
  132.  
  133.     put ""
  134.     put "eigenvalues for array 4:"
  135.     put ""
  136.  
  137.     eigenvv( A, iwork, z, wr, wi, cnt )
  138.  
  139.     for i := 0...DIM-1 do
  140.         put complex_str( wr[i], wi[i] ), "   it = ", cnt[i]
  141.     end for
  142.  
  143. end program
  144.  
  145. %
  146. % compute eigenvalues and eigenvectors of a matrix
  147. %
  148. procedure eigenvv( var a : rmatrix, 
  149.                    var iwork : ivector, 
  150.                    var z : rmatrix, 
  151.                    var wr, wi : rvector, 
  152.                    var cnt : ivector )
  153.  
  154.     var ierr : int
  155.  
  156.     elmhes( a, iwork )
  157.  
  158.     eltran( a, iwork, z )
  159.  
  160.     hqr2( a, wr, wi, z, ierr, cnt )
  161.  
  162.     if ierr ~= 0 then
  163.  
  164.         put "ierr = ", ierr
  165.  
  166.     end if
  167.  
  168. end procedure
  169.  
  170. %
  171. % convert a general matrix to Hessenberg form
  172. %
  173. procedure elmhes( var a : rmatrix, var intar : ivector )
  174.  
  175.     var i, j, m, la, kp1, mm1, mp1 : int
  176.     var x, y : real
  177.  
  178.     la := DIM - 2
  179.     kp1 := 1
  180.  
  181.     if la < kp1 then
  182.         return
  183.     end if
  184.  
  185.     for m := kp1...la do
  186.  
  187.         mm1 := m - 1
  188.         x := 0.0
  189.         i := m
  190.  
  191.         for j := m...DIM-1 do
  192.  
  193.             if rabs( a[j,mm1] ) > rabs( x ) then
  194.  
  195.                 x := a[j,mm1]
  196.                 i := j
  197.  
  198.             end if
  199.  
  200.         end for
  201.  
  202.         intar[m] := i
  203.  
  204.         if i ~= m then
  205.  
  206.             for j := mm1...DIM-1 do
  207.  
  208.                 y := a[i,j]
  209.                 a[i,j] := a[m,j]
  210.                 a[m,j] := y
  211.  
  212.             end for
  213.  
  214.             for j := 0...DIM-1 do
  215.  
  216.                 y := a[j,i]
  217.                 a[j,i] := a[j,m]
  218.                 a[j,m] := y
  219.  
  220.             end for
  221.  
  222.         end if
  223.  
  224.         if x ~= 0.0 then
  225.  
  226.             mp1 := m + 1
  227.  
  228.             for i := mp1...DIM-1 do
  229.  
  230.                 y := a[i,mm1]
  231.                 
  232.                 if y ~= 0.0 then
  233.  
  234.                     y := y / x
  235.                     a[i,mm1] := y
  236.  
  237.                     for j := m...DIM-1 do
  238.  
  239.                         a[i,j] := a[i,j] - y * a[m,j]
  240.                     
  241.                     end for
  242.  
  243.                     for j := 0...DIM-1 do
  244.  
  245.                         a[j,m] := a[j,m] + y * a[j,i]
  246.  
  247.                     end for
  248.  
  249.                 end if
  250.  
  251.             end for
  252.  
  253.         end if
  254.  
  255.     end for
  256.  
  257. end procedure
  258.    
  259. %
  260. % create transform matrix using results of elmhes
  261. %
  262. procedure eltran( var a : rmatrix, 
  263.                   var intar : ivector, 
  264.                   var z : rmatrix )
  265.  
  266.     var i, j, k, l, kl, mn, mp, mp1 : int
  267.  
  268.     % set z to identity matrix
  269.     for i := 0...DIM-1 do
  270.  
  271.         for j := 0...DIM-1 do
  272.  
  273.             z[i,j] := 0.0
  274.  
  275.         end for
  276.  
  277.         z[i,i] := 1.0
  278.  
  279.     end for
  280.  
  281.     kl := DIM - 2
  282.  
  283.     if kl < 1 then
  284.  
  285.         return
  286.  
  287.     end if
  288.  
  289.     for decreasing i := DIM-2...1 do
  290.  
  291.         j := intar[i]
  292.         for k := i+1...DIM-1 do
  293.  
  294.             z[k,i] := a[k,i-1]
  295.  
  296.         end for
  297.  
  298.         if i ~= j then
  299.  
  300.             for k := i...DIM-1 do
  301.  
  302.                 z[i,k] := z[j,k]
  303.                 z[j,k] := 0.0
  304.  
  305.             end for
  306.  
  307.             z[j,i] := 1.0
  308.  
  309.         end if
  310.  
  311.     end for
  312.  
  313. end procedure
  314.  
  315. %
  316. % apply QR method to Hessenberg matrix
  317. %
  318. procedure hqr2( var h : rmatrix, 
  319.                 var wr, wi : rvector, 
  320.                 var vecs : rmatrix, 
  321.                 var ierr : int, 
  322.                 var cnt : ivector )
  323.  
  324.     var i, j, k, l, m, ii, jj, kk, mm, na, nm, nn, its, mpr, enm2, en : int
  325.     var it1, it2, itlim : int
  326.     var p, q, r, s, t, w, x, y, z, ra, sa, vi, vr, norm : real
  327.     var xr, yr, xi, yi, zr, zi : real
  328.     var notlast : boolean
  329.     label nextw:
  330.     label nextit:
  331.     label break1:
  332.  
  333.     it1 := 10
  334.     it2 := 20
  335.     itlim := 30
  336.     ierr := 0
  337.     
  338.     norm := 0.0
  339.     k := 0
  340.  
  341.     for i := 0...DIM-1 do
  342.  
  343.         for j := k...DIM-1 do
  344.  
  345.             norm := norm + rabs( h[i,j] )
  346.  
  347.         end for
  348.  
  349.         k := i
  350.         cnt[i] := 0
  351.  
  352.     end for
  353.  
  354.     assert norm > 0.0
  355.  
  356.     en := DIM-1
  357.     t := 0.0
  358.  
  359. nextw:
  360.  
  361.     watch( en )
  362.  
  363.     if en >= 0 then
  364.  
  365.         its := 0
  366.         na := en - 1
  367.         enm2 := na - 1
  368.  
  369. nextit:
  370.  
  371.         for decreasing l := en...1 do
  372.  
  373.             s := rabs( h[l-1,l-1] ) + rabs( h[l,l] )
  374.  
  375.             if s = 0.0 then
  376.  
  377.                 s := norm
  378.  
  379.             end if
  380.  
  381.             if rabs( h[l,l-1] ) <= s * macheps then
  382.  
  383.                 goto break1
  384.  
  385.             end if
  386.  
  387.         end for
  388.  
  389.         l := 0
  390.  
  391. break1:
  392.  
  393.         x := h[en,en]
  394.  
  395.         watch( l )
  396.  
  397.         if l = en then      % found single a root
  398.  
  399.             wr[en] := x + t
  400.             wi[en] := 0.0
  401.             h[en,en] := x + t
  402.             cnt[en] := its
  403.             en := na
  404.             goto nextw
  405.  
  406.         end if
  407.  
  408.         y := h[na,na]
  409.         w := h[en,na] * h[na,en]
  410.  
  411.         if l ~= na then
  412.  
  413.             if its = itlim then
  414.  
  415.                 cnt[en] := 31
  416.                 ierr := en
  417.                 return
  418.  
  419.             end if
  420.  
  421.             if its = it1 or its = it2 then
  422.  
  423.                 t := t + x
  424.                 for i := 0...en do
  425.  
  426.                     h[i,i] := h[i,i] - x
  427.  
  428.                 end for
  429.  
  430.                 s := rabs( h[en,na] ) + rabs( h[na,en-2] )
  431.                 x := 0.75 * s
  432.                 y := x
  433.                 w := -0.4375 * s * s
  434.  
  435.             end if
  436.  
  437.             incr its
  438.  
  439.             for decreasing m := en-2...l do
  440.  
  441.                 mm := m + 1
  442.                 z := h[m,m]
  443.                 r := x - z
  444.                 s := y - z
  445.                 p := ( r * s - w ) / ( h[mm,m] ) + h[m,mm]
  446.                 q := h[mm,mm] - z - r - s
  447.                 r := h[m+2,mm]
  448.                 s := rabs( p ) + rabs( q ) + rabs( r )
  449.                 p := p / s
  450.                 q := q / s
  451.                 r := r / s
  452.  
  453.                 if m = l then
  454.  
  455.                     exit
  456.  
  457.                 end if
  458.  
  459.                 if rabs( h[m,m-1] ) * ( rabs( q ) + rabs( r ) ) <=
  460.                    macheps * rabs( p ) * 
  461.                    ( rabs( h[m-1,m-1] ) + rabs( z ) + rabs( h[mm,mm] ) )
  462.                 then
  463.  
  464.                     exit
  465.  
  466.                 end if
  467.  
  468.             end for
  469.  
  470.             for i := m+2...en do
  471.  
  472.                 h[i,i-2] := 0.0
  473.  
  474.             end for
  475.  
  476.             for i := m+3...en do
  477.  
  478.                 h[i,i-3] := 0.0
  479.  
  480.             end for
  481.  
  482.             for k := m...na do
  483.  
  484.                 notlast := k ~= na
  485.  
  486.                 if k ~= m then
  487.  
  488.                     p := h[k,k-1]
  489.                     q := h[k+1,k-1]
  490.                     
  491.                     if notlast then
  492.  
  493.                         r := h[k+2,k-1]
  494.  
  495.                     else
  496.  
  497.                         r := 0.0
  498.  
  499.                     end if
  500.  
  501.                     x := rabs( r ) + rabs( q ) + rabs( p )
  502.  
  503.                     if x = 0.0 then
  504.  
  505.                         exit
  506.  
  507.                     end if
  508.  
  509.                     p := p / x
  510.                     q := q / x
  511.                     r := r / x
  512.  
  513.                 end if
  514.  
  515.                 s := sqrt( p * p + q * q + r * r )
  516.  
  517.                 if p < 0.0 then
  518.  
  519.                     s := -s
  520.  
  521.                 end if
  522.  
  523.                 if k ~= m then
  524.  
  525.                     h[k,k-1] := -s * x
  526.  
  527.                 elsif l ~= m then
  528.  
  529.                     h[k,k-1] := -h[k,k-1]
  530.  
  531.                 end if
  532.  
  533.                 p := p + s
  534.  
  535.                 x := p / s
  536.                 y := q / s
  537.                 z := r / s
  538.                 
  539.                 q := q / p
  540.                 r := r / p
  541.  
  542.                 % row modification
  543.  
  544.                 for j := k...DIM-1 do
  545.  
  546.                     p := h[k,j] + q * h[k+1,j]
  547.  
  548.                     if notlast then
  549.  
  550.                         p := p + r * h[k+2,j]
  551.                         h[k+2,j] := h[k+2,j] - p * z
  552.  
  553.                     end if
  554.  
  555.                     h[k+1,j] := h[k+1,j] - p * y
  556.                     h[k,j] := h[k,j] - p * x
  557.  
  558.                 end for
  559.  
  560.                 % column modification
  561.  
  562.                 if k + 3 < en then
  563.  
  564.                     j := k + 3
  565.  
  566.                 else
  567.  
  568.                     j := en
  569.  
  570.                 end if
  571.  
  572.                 for i := 0...j do
  573.  
  574.                     p := x * h[i,k] + y * h[i,k+1]
  575.  
  576.                     if notlast then
  577.  
  578.                         p := p + z * h[i,k+2]
  579.                         h[i,k+2] := h[i,k+2] - p * r
  580.  
  581.                     end if
  582.  
  583.                     h[i,k+1] := h[i,k+1] - p * q
  584.                     h[i,k] := h[i,k] - p
  585.  
  586.                 end for
  587.  
  588.                 % accumulate
  589.  
  590.                 for i := 0...DIM-1 do
  591.  
  592.                     p := x * vecs[i,k] + y * vecs[i,k+1]
  593.                     
  594.                     if notlast then
  595.  
  596.                         p := p + z * vecs[i,k+2]
  597.                         vecs[i,k+2] := vecs[i,k+2] - p * r
  598.  
  599.                     end if
  600.  
  601.                     vecs[i,k+1] := vecs[i,k+1] - p * q
  602.                     vecs[i,k] := vecs[i,k] - p
  603.  
  604.                 end for
  605.  
  606.             end for % k
  607.  
  608.             goto nextit
  609.  
  610.         end if % l ~= na
  611.  
  612.         % else two roots were found
  613.  
  614.         p := 0.5 * ( y - x )
  615.         q := p * p + w
  616.         z := sqrt( rabs( q ) )
  617.         x := x + t
  618.         h[en,en] := x
  619.         h[na,na] := y + t
  620.         cnt[en] := -its
  621.         cnt[na] :=  its
  622.  
  623.         if q >= 0.0 then    % real pair of roots
  624.  
  625.             if p < 0.0 then
  626.  
  627.                 z := p - z
  628.  
  629.             else
  630.  
  631.                 z := p + z
  632.  
  633.             end if
  634.  
  635.             wr[na] := x + z
  636.             wr[en] := wr[na]
  637.  
  638.             if z ~= 0.0 then
  639.  
  640.                 s := x - w / z
  641.                 wr[en] := s
  642.  
  643.             end if
  644.  
  645.             wi[en] := 0.0
  646.             wi[na] := 0.0
  647.  
  648.             x := h[en,na]
  649.             s := rabs( x ) + rabs( z )
  650.             p := x / s
  651.             q := z / s
  652.             r := sqrt( p * p + q * q )
  653.             p := p / r
  654.             q := q / r
  655.  
  656.             for j := na...DIM-1 do
  657.  
  658.                 z := h[na,j]
  659.                 h[na,j] := q * z + p * h[en,j]
  660.                 h[en,j] := q * h[en,j] - p * z
  661.  
  662.             end for
  663.  
  664.             for i := 0...en do
  665.  
  666.                 z := h[i,na]
  667.                 h[i,na] := q * z + p * h[i,en]
  668.                 h[i,en] := q * h[i,en] - p * z
  669.                 
  670.             end for
  671.  
  672.             for i := 0...DIM-1 do
  673.  
  674.                 z := vecs[i,na]
  675.                 vecs[i,na] := q * z + p * vecs[i,en]
  676.                 vecs[i,en] := q * vecs[i,en] - p * z
  677.             
  678.             end for
  679.  
  680.         else                % complex pair of roots
  681.  
  682.             wr[en] := x + p
  683.             wr[na] := x + p
  684.             wi[na] := z
  685.             wi[en] := -z
  686.  
  687.         end if
  688.  
  689.         decr en
  690.         decr en
  691.  
  692.         goto nextw
  693.  
  694.     end if % en >= 0
  695.  
  696.     % all eigenvalues have been found
  697.  
  698. end procedure
  699.  
  700.  
  701. %=====================================================================
  702. %                     some utility functions
  703. %=====================================================================
  704.  
  705. %
  706. % return the absolute value of a real number
  707. %
  708. function rabs( x : real ) : real
  709.  
  710.     var rv : real
  711.  
  712.     if x < 0.0 then
  713.  
  714.         rv := -x
  715.  
  716.     else
  717.  
  718.         rv := x
  719.  
  720.     end if
  721.  
  722.     return rv
  723.  
  724. end function
  725.  
  726. %
  727. % print a real 2D matrix
  728. %
  729. procedure print_matrix( var a : rmatrix )
  730.  
  731.     var i, j, btm, top, count : int := 0
  732.  
  733.     put ""
  734.     
  735.     loop
  736.  
  737.         exit when btm >= DIM
  738.  
  739.         if DIM > btm + 8 then
  740.  
  741.             top := btm + 8
  742.  
  743.         else
  744.  
  745.             top := DIM
  746.  
  747.         end if
  748.  
  749.  
  750.         put "matrix columns ", btm, " to ", top - 1
  751.  
  752.         for j := 0...DIM-1 do
  753.  
  754.             for i := btm...top-1 do
  755.  
  756.                 put a[j,i]:15:6:3...
  757.  
  758.             end for
  759.  
  760.             put ""
  761.  
  762.         end for
  763.  
  764.         btm := btm + 8
  765.  
  766.     end loop
  767.  
  768. end procedure
  769.  
  770. %
  771. % convert complex number to a string
  772. %
  773. function complex_str( re, im : real ) : string
  774.  
  775.     var s : string
  776.  
  777.     if im >= 0.0 then
  778.  
  779.         s := erealstr(re,15,6,3) & " +" & erealstr(im,13,6,3) & "j"
  780.  
  781.     else
  782.  
  783.         im := -im
  784.         s := erealstr(re,15,6,3) & " -" & erealstr(im,13,6,3) & "j"
  785.  
  786.     end if
  787.  
  788.     return s
  789.  
  790. end function